#####################################################################
#Program for fit our generalization on Hausman`s approach on \pi_1=probit and \pi_0=0.(Our METHOD 2)
#####################################################################

misclass_phi1=function (form,  bmat = 0, print.summary = TRUE, 
                        data = NULL) # Here form must be y_both~x+z1+z2
{
  probit <- glm(form, family = binomial(link = "probit"), data = data) 
  #cat("Standard Probit Estimates", "\n")
  #print(summary(probit))
  ysum <- model.frame(form, data = data)[, 1]
  xzmat <- model.matrix(form, data = data) # Design matrix (1,x,z1,z2)
  
  probitX <- glm(ysum~ xzmat[,2], family = binomial(link = "probit"), data = data)  
  xmat<- model.matrix(ysum~ xzmat[,2]) # Design matrix (1,x)
  if (identical(bmat,0)){
    nbxzmat <- length(probit$coef)
    bxmat <- probitX$coef
    bxzmat <- probit$coef
    bmat<-c(bxmat,bxzmat)
  }
  nk = length(probit$coef)+ length(probitX$coef)
  bstart <-  bmat
  names(bstart)<-c("beta0","beta1","eta0","eta1","eta2","eta3")
  
  logl <- function(b) {
    eysum<-(1-pnorm(xzmat%*% b[3:nk]))*pnorm(xmat%*%b[1:2]) 
    out<--sum(ifelse(ysum==1,log(eysum),log(1-eysum)))
    return(out)
  }
  # nlfit<- optim( bstart,logl,method="L-BFGS-B", hessian = TRUE)
  #  nlfit<- nlm( logl, starts, hessian = TRUE, iterlim = 1000)
  
  nlfit<-optim(bstart,logl, gr=NULL, method="BFGS",control=list(maxit=500), hessian=TRUE) # the function probit.lf is defined above
  # For default maxit=100 and the iteration limit may need to be larger to get convergence. 
  vmat <- solve(nlfit$hessian)
  stderr <- sqrt(abs(diag(vmat)))
  tval <- nlfit$par/stderr
  pval <- 2 * (1 - pnorm(abs(tval)))
  
  cname <- c("Estimate", "Std. Error", "z-value", "Pr(>|z|)")
  if (print.summary == TRUE) {
    outmat <- cbind(nlfit$par, stderr, tval, pval)
    #rownames(outmat) <- names(probit$coef)
    colnames(outmat) <- cname
    print(round(outmat, 5))
    cat(" ", "\n")
  }
  
  out1 <- list(nlfit$par, stderr, vmat, nlfit$convergence, 
               -nlfit$value)
  names(out1) <- c( "estimate", "stderr", "vmat","convergence", "minimum")
  return(out1)  
}
